/* Copyright (c) 2021, National University of Defense Technology. All rights reserved.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 * 
 *     http://www.apache.org/licenses/LICENSE-2.0
 * 
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include <mpi.h>

namespace YHAMG
{

#ifndef VECTOR_H
#define VECTOR_H

struct Vector
{
	int ref;
	int size;
	double* values;

	Vector();
	Vector(int n, double* values, int ref);
	Vector(const Vector& x);
	Vector(Vector&& x);
	~Vector();
	Vector& operator=(double a);
	Vector& operator=(const Vector& x);
	Vector& operator=(Vector&& x);

	void Free();
	void Resize(int n);
	void Fill(double a) const;
	void FillRandom() const;
	void Copy(const Vector& x) const;
	void Scale(double a) const;
	void Shift(double a) const;
	void AddScaled(double a, const Vector& x) const;
	void Add2Scaled(double a, const Vector& x, double b, const Vector& y) const;
	void Reference(const Vector& x);
};
#endif

#ifndef MULTIVECTOR_H
#define MULTIVECTOR_H

struct MultiVector
{
	int ref;
	int num_vectors;
	int size;
	double* values;

	MultiVector();
	MultiVector(int m, int n, double* values, int ref);
	MultiVector(const MultiVector& X);
	MultiVector(MultiVector&& X);
	~MultiVector();
	MultiVector& operator=(const MultiVector& X);
	MultiVector& operator=(MultiVector&& X);

	void Free();
	void Allocate(int m, int n);
	void Reference(const MultiVector& X);
	const Vector operator()(int j) const;
};

#endif

#ifndef PAR_VECTOR_H
#define PAR_VECTOR_H

struct Par_Vector
{
	MPI_Comm comm;
	Vector Loc;

	Par_Vector();
	Par_Vector(MPI_Comm comm);
	Par_Vector(MPI_Comm comm, int Loc_size, double* Loc_values, int Loc_ref);
	Par_Vector(const Par_Vector& x);
	Par_Vector(Par_Vector&& x);
	Par_Vector& operator=(double a);
	Par_Vector& operator=(const Par_Vector& x);
	Par_Vector& operator=(Par_Vector&& x);

	void Free();
	void Resize(int n);
	void Reference(const Par_Vector& x);
	void Fill(double a) const;
	void FillRandom() const;
	void Copy(const Par_Vector& x) const;
	void Scale(double a) const;
	void Shift(double a) const;
	void AddScaled(double a, const Par_Vector& x) const;
	void Add2Scaled(double a, const Par_Vector& x, double b, const Par_Vector& y) const;
};

#endif

#ifndef PAR_MULTIVECTOR_H
#define PAR_MULTIVECTOR_H

struct Par_MultiVector
{
	MPI_Comm comm;
	int num_vectors;
	MultiVector Loc;

	Par_MultiVector();
	Par_MultiVector(MPI_Comm comm);
	Par_MultiVector(MPI_Comm comm, int num_vectors, int Loc_size, double* Loc_values, int Loc_ref);
	Par_MultiVector(const Par_MultiVector& X);
	Par_MultiVector(Par_MultiVector&& X);
	Par_MultiVector& operator=(const Par_MultiVector& X);
	Par_MultiVector& operator=(Par_MultiVector&& X);

	void Free();
	void Allocate(int m, int n);
	void Reference(const Par_MultiVector& X);
	const Par_Vector operator()(int j) const;
};

#endif

#ifndef PARVEC_DOT_H
#define PARVEC_DOT_H

double ParVec_Dot(const Par_Vector& x, const Par_Vector& y);

#endif

#ifndef PARVEC_AXPBY_H
#define PARVEC_AXPBY_H

void ParVec_AXPBY(double alpha, const Par_Vector& x, double beta, const Par_Vector& y, const Par_Vector& w);

#endif

#ifndef PARVEC_AXPBYPCZ_H
#define PARVEC_AXPBYPCZ_H

void ParVec_AXPBYPCZ(double alpha, const Par_Vector& x, double beta, const Par_Vector& y, double gamma, const Par_Vector& z, const Par_Vector& w);

#endif

#ifndef OPERATOR_H
#define OPERATOR_H

class Operator
{
public:
	virtual ~Operator();
	virtual int InSize() const = 0;
	virtual int OutSize() const = 0;
	virtual void Apply(const Vector& x, const Vector& y) const = 0;
	virtual void Apply(const MultiVector& X, const MultiVector& Y) const;
};

#endif

#ifndef PAR_OPERATOR_H
#define PAR_OPERATOR_H

class Par_Operator
{
public:
	MPI_Comm comm;

	Par_Operator();
	Par_Operator(MPI_Comm comm);
	virtual ~Par_Operator();
	virtual int InSize() const = 0;
	virtual int OutSize() const = 0;
	virtual void Apply(const Par_Vector& x, const Par_Vector& y) const = 0;
	virtual void Apply(const Par_MultiVector& X, const Par_MultiVector& Y) const;
};

#endif

#ifndef PAR_ZERO_OPERATOR_H
#define PAR_ZERO_OPERATOR_H

class Par_Zero_Operator : public Par_Operator
{
private:
	int size[2];

public:
	Par_Zero_Operator(MPI_Comm comm, int n);
	Par_Zero_Operator(MPI_Comm comm, int n, int m);	
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& x, const Par_Vector& y) const;
};

#endif

#ifndef PAR_IDENTITY_OPERATOR_H
#define PAR_IDENTITY_OPERATOR_H

class Par_Identity_Operator : public Par_Operator
{
private:
	int size[2];

public:
	Par_Identity_Operator(MPI_Comm comm, int n);
	Par_Identity_Operator(MPI_Comm comm, int n, int m);	
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& x, const Par_Vector& y) const;
};

#endif

#ifndef PAR_COMPOUND_OPERATOR_H
#define PAR_COMPOUND_OPERATOR_H

class Par_Compound_Operator : public Par_Operator
{
private:
	const Par_Operator* A;
	const Par_Operator* B;
	const Par_Operator* C;
	double alpha;
	double beta;
	double gamma;

public:	
	Par_Compound_Operator(double alpha, const Par_Operator& A, double beta);	
	Par_Compound_Operator(double alpha, const Par_Operator& A, double beta, const Par_Operator& B, double gamma);
	Par_Compound_Operator(double alpha, const Par_Operator& A, const Par_Operator& B, double gamma);
	Par_Compound_Operator(double alpha, const Par_Operator& A, const Par_Operator& B, double beta, const Par_Operator& C, double gamma);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& x, const Par_Vector& y) const;
	void Apply(const Par_MultiVector& X, const Par_MultiVector& Y) const;
};

#endif

#ifndef PAR_CHEBYSHEV_POLYNOMIAL_H
#define PAR_CHEBYSHEV_POLYNOMIAL_H

class Par_Chebyshev_Polynomial : public Par_Operator
{
private:
	const Par_Operator* A;
	double alpha;
	double beta;
	int degree;

public:
	Par_Chebyshev_Polynomial(const Par_Operator& A, double alpha, double beta, int degree);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& x, const Par_Vector& y) const;
	void Apply(const Par_MultiVector& X, const Par_MultiVector& Y) const;
};

#endif

#ifndef PAR_SOLVER_H
#define PAR_SOLVER_H

class Par_Solver
{
public:
	virtual ~Par_Solver();
	virtual void operator()(const Par_Vector& b, const Par_Vector& x, int& iter, double& relres) const = 0;
};

#endif

#ifndef PAR_CG_SOLVER_H
#define PAR_CG_SOLVER_H

class Par_CG_Solver : public Par_Solver
{
private:
	const Par_Operator* A;
	const Par_Operator* P;

public:
	int    MaxIters;
	double Tolerance;
	int    PrintStats;

	Par_CG_Solver(const Par_Operator& A);
	Par_CG_Solver(const Par_Operator& A, const Par_Operator& P);
	Par_CG_Solver(const Par_Operator& A, int max_Iters, double tolerance, int print_stats);
	Par_CG_Solver(const Par_Operator& A, const Par_Operator& P, int max_Iters, double tolerance, int print_stats);
	void operator()(const Par_Vector& b, const Par_Vector& x, int& iter, double& relres) const;
};


#endif

#ifndef PAR_GMRES_SOLVER_H
#define PAR_GMRES_SOLVER_H

class Par_GMRES_Solver : public Par_Solver
{
private:
	const Par_Operator* A;
	const Par_Operator* P;

public:
	int    Restart;
	int    MaxIters;
	double Tolerance;
	int    PrintStats;

	Par_GMRES_Solver(const Par_Operator& A);
	Par_GMRES_Solver(const Par_Operator& A, const Par_Operator& P);
	Par_GMRES_Solver(const Par_Operator& A, int restart, int max_Iters, double tolerance, int print_stats);
	Par_GMRES_Solver(const Par_Operator& A, const Par_Operator& P, int restart, int max_Iters, double tolerance, int print_stats);
	void operator()(const Par_Vector& b, const Par_Vector& x, int& iter, double& relres) const;
};

#endif

#ifndef PAR_BICGSTAB_SOLVER_H
#define PAR_BICGSTAB_SOLVER_H

class Par_BiCGStab_Solver : public Par_Solver
{
private:
	const Par_Operator* A;
	const Par_Operator* P;

public:
	int    MaxIters;
	double Tolerance;
	int    PrintStats;

	Par_BiCGStab_Solver(const Par_Operator& A);
	Par_BiCGStab_Solver(const Par_Operator& A, const Par_Operator& P);
	Par_BiCGStab_Solver(const Par_Operator& A, int max_Iters, double tolerance, int print_stats);
	Par_BiCGStab_Solver(const Par_Operator& A, const Par_Operator& P, int max_Iters, double tolerance, int print_stats);
	void operator()(const Par_Vector& b, const Par_Vector& x, int& iter, double& relres) const;
};

#endif

#ifndef PAR_PIPECG_SOLVER_H
#define PAR_PIPECG_SOLVER_H

class Par_PipeCG_Solver : public Par_Solver
{
private:
	const Par_Operator* A;
	const Par_Operator* P;

public:
	int    MaxIters;
	double Tolerance;
	int    PrintStats;

	Par_PipeCG_Solver(const Par_Operator& A);
	Par_PipeCG_Solver(const Par_Operator& A, const Par_Operator& P);
	Par_PipeCG_Solver(const Par_Operator& A, int max_Iters, double tolerance, int print_stats);
	Par_PipeCG_Solver(const Par_Operator& A, const Par_Operator& P, int max_Iters, double tolerance, int print_stats);
	void operator()(const Par_Vector& b, const Par_Vector& x, int& iter, double& relres) const;
};

#endif

#ifndef PAR_PIPEGMRES_SOLVER_H
#define PAR_PIPEGMRES_SOLVER_H

class Par_PipeGMRES_Solver : public Par_Solver
{
private:
	const Par_Operator* A;
	const Par_Operator* P;

public:
	int    Restart;
	int    MaxIters;
	double Shift;
	double Tolerance;
	int    PrintStats;

	Par_PipeGMRES_Solver(const Par_Operator& A);
	Par_PipeGMRES_Solver(const Par_Operator& A, const Par_Operator& P);
	Par_PipeGMRES_Solver(const Par_Operator& A, int restart, double shift, int max_Iters, double tolerance, int print_stats);
	Par_PipeGMRES_Solver(const Par_Operator& A, const Par_Operator& P, int restart, double shift, int max_Iters, double tolerance, int print_stats);
	void operator()(const Par_Vector& b, const Par_Vector& x, int& iter, double& relres) const;
};

#endif

#ifndef PAR_PIPEBICGSTAB_SOLVER_H
#define PAR_PIPEBICGSTAB_SOLVER_H

class Par_PipeBiCGStab_Solver : public Par_Solver
{
private:
	const Par_Operator* A;
	const Par_Operator* P;

public:
	int    MaxIters;
	double Tolerance;
	int    PrintStats;

	Par_PipeBiCGStab_Solver(const Par_Operator& A);
	Par_PipeBiCGStab_Solver(const Par_Operator& A, const Par_Operator& P);
	Par_PipeBiCGStab_Solver(const Par_Operator& A, int max_Iters, double tolerance, int print_stats);
	Par_PipeBiCGStab_Solver(const Par_Operator& A, const Par_Operator& P, int max_Iters, double tolerance, int print_stats);
	void operator()(const Par_Vector& b, const Par_Vector& x, int& iter, double& relres) const;
};

#endif

#ifndef PAR_EIGENSOLVER_H
#define PAR_EIGENSOLVER_H

class Par_EigenSolver
{
public:
	virtual ~Par_EigenSolver();
	virtual void operator()(const Par_MultiVector& X, const Vector& lamdba, int& iter, const Vector& res) const = 0;
};

#endif

#ifndef PAR_LOBPCG_EIGENSOLVER_H
#define PAR_LOBPCG_EIGENSOLVER_H

class Par_LOBPCG_EigenSolver : public Par_EigenSolver
{
private:
	const Par_Operator* A;
	const Par_Operator* B;
	const Par_Operator* T;
	const Par_MultiVector* Y;
	
public:
	int    MaxIters;
	double Tolerance;
	int    PrintStats;

	Par_LOBPCG_EigenSolver(const Par_Operator& A);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_MultiVector& Y);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_Operator& T);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_Operator& T, const Par_MultiVector& Y);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_Operator& B, const Par_Operator& T);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_Operator& B, const Par_Operator& T, const Par_MultiVector& Y);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, int max_Iters, double tolerance, int print_stats);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_MultiVector& Y, int max_Iters, double tolerance, int print_stats);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_Operator& T, int max_Iters, double tolerance, int print_stats);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_Operator& T, const Par_MultiVector& Y, int max_Iters, double tolerance, int print_stats);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_Operator& B, const Par_Operator& T, int max_Iters, double tolerance, int print_stats);
	Par_LOBPCG_EigenSolver(const Par_Operator& A, const Par_Operator& B, const Par_Operator& T, const Par_MultiVector& Y, int max_Iters, double tolerance, int print_stats);
	void operator()(const Par_MultiVector& X, const Vector& lamdba, int& iter, const Vector& res) const;
};

#endif

#ifndef COO_MATRIX_H
#define COO_MATRIX_H

struct COO_Matrix : public Operator
{
	int ref;
	int size[2];
	int nnz;
	
	int* rowind;
	int* colind;
	double* values;

	COO_Matrix();
	COO_Matrix(int n, int m, int nnz, int* rowind, int* colind, double* values, int ref);
	COO_Matrix(const COO_Matrix& A);
	COO_Matrix(COO_Matrix&& A);
	~COO_Matrix();
	COO_Matrix& operator=(const COO_Matrix& A);
	COO_Matrix& operator=(COO_Matrix&& A);
	
	void Free();
	void Reference(const COO_Matrix& A);
	int InSize() const;
	int OutSize() const;
	void Apply(const Vector& x, const Vector& y) const;
};

#endif

#ifndef CSR_MATRIX_H
#define CSR_MATRIX_H

struct CSR_Matrix : public Operator
{
	int ref;
	int size[2];

	int* rowptr;
	int* colind;
	double* values;

	CSR_Matrix();
	CSR_Matrix(int n, int m, int* rowptr, int* colind, double* values, int ref);
	CSR_Matrix(const CSR_Matrix& A);
	CSR_Matrix(CSR_Matrix&& A);
	CSR_Matrix(const COO_Matrix& A);
	~CSR_Matrix();
	CSR_Matrix& operator=(const CSR_Matrix& A);
	CSR_Matrix& operator=(CSR_Matrix&& A);
	CSR_Matrix& operator=(const COO_Matrix& A);

	void Free();
	void Reference(const CSR_Matrix& A);
	int InSize() const;
	int OutSize() const;
	void Apply(const Vector& x, const Vector& y) const;
	void Apply(const MultiVector& X, const MultiVector& Y) const;
};


#endif

#ifndef PARCSR_MATRIX_H
#define PARCSR_MATRIX_H

struct ParCSR_Matrix : public Par_Operator
{
	CSR_Matrix LocLoc;
	CSR_Matrix LocExt;

	int num_neighbors;
	int* neighbor_ranks;

	int* Ext_ptr;
	int* Ext_map;

	int* send_ptr;
	int* send_list;
	double* send_buffer;

	Vector x_Ext;

	ParCSR_Matrix();
	ParCSR_Matrix(MPI_Comm comm, int LocLoc_rows, int LocLoc_cols, int LocExt_cols, 
	int* LocLoc_rowptr, int* LocLoc_colind, double* LocLoc_values, int LocLoc_ref, 
	int* LocExt_rowptr, int* LocExt_colind, double* LocExt_values, int LocExt_ref,
	int num_neighbors, int* neighbor_ranks, int* Ext_ptr, int* Ext_map);
	ParCSR_Matrix(const ParCSR_Matrix& A);
	ParCSR_Matrix(ParCSR_Matrix&& A);
	~ParCSR_Matrix();
	ParCSR_Matrix& operator=(const ParCSR_Matrix& A);
	ParCSR_Matrix& operator=(ParCSR_Matrix&& A);
	
	void Free();
	void Reference(const ParCSR_Matrix& A);
	void ExchangeHalo(const Vector& x_Loc, const Vector& x_Ext) const;
	void SetupHalo();
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& x, const Par_Vector& y) const;
	void Apply(const Par_MultiVector& X, const Par_MultiVector& Y) const;
};

#endif

#ifndef PARCSR_MATVEC_H
#define PARCSR_MATVEC_H

void ParCSR_MatVec(double alpha, const ParCSR_Matrix& A, const Par_Vector& x, double beta, const Par_Vector& y, const Par_Vector& w);
void ParCSR_MatVec(const ParCSR_Matrix& A, const Par_Vector& x, const Par_Vector& y);

#endif

#ifndef PARCSR_MATMUL_H
#define PARCSR_MATMUL_H

void ParCSR_MatMul(const ParCSR_Matrix& A, const ParCSR_Matrix& B, ParCSR_Matrix& C);

#endif

#ifndef PARCSR_TRANSPOSE_H
#define PARCSR_TRANSPOSE_H

void ParCSR_Transpose(const ParCSR_Matrix& A, ParCSR_Matrix& B);

#endif

#ifndef PARCSR_GENERATOR_H
#define PARCSR_GENERATOR_H

class ParCSR_Generator
{
public:
	MPI_Comm comm;

	ParCSR_Generator();
	ParCSR_Generator(MPI_Comm comm);
	virtual ~ParCSR_Generator();
	virtual void operator()(ParCSR_Matrix& A) const = 0;
};

#endif

#ifndef PARCSR_LAPLACIAN_Generator_H
#define PARCSR_LAPLACIAN_Generator_H

class ParCSR_Laplacian_Generator : public ParCSR_Generator
{
public:
	int nx;
	int ny;
	int nz;

	int Px;
	int Py;
	int Pz;

	double cx;
	double cy;
	double cz;

	ParCSR_Laplacian_Generator();
	ParCSR_Laplacian_Generator(MPI_Comm comm);
	ParCSR_Laplacian_Generator(MPI_Comm comm, int nx, int ny, int nz, int Px, int Py, int Pz, double cx, double cy, double cz);
	void operator()(ParCSR_Matrix& A) const;
};

#endif

#ifndef PARCSR_LAPLACIAN27PT_Generator_H
#define PARCSR_LAPLACIAN27PT_Generator_H

class ParCSR_Laplacian27pt_Generator : public ParCSR_Generator
{
public:
	int nx;
	int ny;
	int nz;

	int Px;
	int Py;
	int Pz;

	ParCSR_Laplacian27pt_Generator();
	ParCSR_Laplacian27pt_Generator(MPI_Comm comm);
	ParCSR_Laplacian27pt_Generator(MPI_Comm comm, int nx, int ny, int nz, int Px, int Py, int Pz);
	void operator()(ParCSR_Matrix& A) const;
};

#endif

#ifndef PARCSR_IJ_GENERATOR_H
#define PARCSR_IJ_GENERATOR_H

class ParCSR_IJ_Generator : public ParCSR_Generator
{
public:
	int Loc_rows;
	const long* col_partitioning;
	const int* LocGlo_rowptr;
	const long* LocGlo_colind;
	const double* LocGlo_values;
	int pattern_symmetry;

	ParCSR_IJ_Generator();
	ParCSR_IJ_Generator(MPI_Comm comm);
	ParCSR_IJ_Generator(MPI_Comm comm, int Loc_rows, const long* col_partitioning, 	const int* LocGlo_rowptr, const long* LocGlo_colind, const double* LocGlo_values, int pattern_symmetry);
	void operator()(ParCSR_Matrix& A) const;
};

#endif

#ifndef PARCSR_PRECONDITIONER_H
#define PARCSR_PRECONDITIONER_H

class ParCSR_Preconditioner : public Par_Operator
{
public:
	virtual ~ParCSR_Preconditioner();
	virtual void Setup(const ParCSR_Matrix& A, int REUSE = 0) = 0;
};

#endif

#ifndef PARCSR_JACOBI_PRECONDITIONER_H
#define PARCSR_JACOBI_PRECONDITIONER_H

class ParCSR_Jacobi_Preconditioner : public ParCSR_Preconditioner
{
private:
	Vector D_Recip;

public:
	void Free();
	void Setup(const ParCSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARCSR_SOR_PRECONDITIONER_H
#define PARCSR_SOR_PRECONDITIONER_H

class ParCSR_SOR_Preconditioner : public ParCSR_Preconditioner
{
private:
	Vector D_Loc;
	Vector D_Recip_Loc;
	CSR_Matrix L_LocLoc;
	CSR_Matrix L_LocExt;
	CSR_Matrix U_LocLoc;
	CSR_Matrix U_LocExt;
	Vector x_Ext;

	bool* interior;

	int label;
	int num_neighbors;
	int* neighbor_labels;
	int* neighbor_ranks;

	int* recv_ptr;
	int* send_ptr;
	int* send_list;
	double* send_buffer;

public:
	int Type;
	double Factor;

	ParCSR_SOR_Preconditioner();
	ParCSR_SOR_Preconditioner(int type, int factor);
	~ParCSR_SOR_Preconditioner();

	void Free();
	void Setup(const ParCSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARCSR_ILU_PRECONDITIONER_H
#define PARCSR_ILU_PRECONDITIONER_H

class ParCSR_ILU_Preconditioner : public ParCSR_Preconditioner
{
private:
	Vector D_Recip_Loc;
	CSR_Matrix L_LocLoc;
	CSR_Matrix L_LocExt;
	CSR_Matrix U_LocLoc;
	CSR_Matrix U_LocExt;
	Vector x_Ext;

	bool* interior;

	int label;
	int num_neighbors;
	int* neighbor_labels;
	int* neighbor_ranks;

	int* recv_ptr;
	int* send_ptr;
	int* send_list;
	double* send_buffer;

public:
	ParCSR_ILU_Preconditioner();
	~ParCSR_ILU_Preconditioner();
	
	void Free();
	void Setup(const ParCSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARCSR_BLOCKJACOBISOR_PRECONDITIONER_H
#define PARCSR_BLOCKJACOBISOR_PRECONDITIONER_H

class ParCSR_BlockJacobiSOR_Preconditioner : public ParCSR_Preconditioner
{
private:
	int num_threads;
	Vector D;
	Vector D_Recip;
	CSR_Matrix L;
	CSR_Matrix U;

public:
	int Type;
	double Factor;

	ParCSR_BlockJacobiSOR_Preconditioner();
	ParCSR_BlockJacobiSOR_Preconditioner(int type, int factor);

	void Free();
	void Setup(const ParCSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARCSR_BLOCKJACOBIILU_PRECONDITIONER_H
#define PARCSR_BLOCKJACOBIILU_PRECONDITIONER_H

class ParCSR_BlockJacobiILU_Preconditioner : public ParCSR_Preconditioner
{
private:
	int num_threads;
	Vector D_Recip;
	CSR_Matrix L;
	CSR_Matrix U;

public:
	ParCSR_BlockJacobiILU_Preconditioner();
	
	void Free();
	void Setup(const ParCSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARCSR_AMG_PRECONDITIONER_H
#define PARCSR_AMG_PRECONDITIONER_H

class ParCSR_AMG_Preconditioner : public ParCSR_Preconditioner
{
private:
	int num_levels;
	ParCSR_Matrix* A;
	ParCSR_Matrix* P;
	ParCSR_Matrix* R;
	class Smoother;
	Smoother* S;

	Par_Vector* r;
	Par_Vector* e;
	Par_Vector* g;

	void Strength(const ParCSR_Matrix& A, bool* LocLoc_s, bool* LocExt_s) const;
	void Coarsening(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map, int AGGRESSIVE = 0) const;
	void Interpolation(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map, ParCSR_Matrix& P, int AGGRESSIVE = 0) const;
	void RAP(const ParCSR_Matrix& R, const ParCSR_Matrix& A, const ParCSR_Matrix& P, ParCSR_Matrix& C) const;
	void Sparsification(const ParCSR_Matrix& A, ParCSR_Matrix& B) const;
	void Smooth(int level, const Par_Vector& b, const Par_Vector& x) const;
	void Restriction(int level, const Par_Vector& b, const Par_Vector& x, const Par_Vector& g) const;
	void Prolongation(int level, const Par_Vector& e, const Par_Vector& x) const;
	void V_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;
	void W_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;
	void F_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;

public:
	int    MaxLevels;
	int    CoarseSize;
	double StrengthThreshold;
	int    AggNumLevels;
	int    CoarsenType;
	int    InterpMinElements;
	int    InterpMaxElements;
	double TruncationFactor;
	double SparsificationThreshold;
	int    CycleType;
	int    PreSweeps;
	int    PostSweeps;
	int    CoarseSweeps;
	int    RelaxationType;
	double RelaxationFactor;
	int    PrintStats;

	ParCSR_AMG_Preconditioner();
	ParCSR_AMG_Preconditioner(
	int    max_levels,
	int    coarse_size,
	double strength_threshold,
	int    agg_num_levels,
	int    coarsen_type,
	int    interp_min_elements,
	int    interp_max_elements,	
	double truncation_factor,
	double sparsification_threshold,
	int    cycle_type,
	int    pre_sweeps,
	int    post_sweeps,
	int    coarse_sweeps,
	int    relaxation_type,
	double relaxation_factor,
	int    print_stats);
	~ParCSR_AMG_Preconditioner();
	
	void Free();
	void Setup(const ParCSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARCSR_ILUAMG_PRECONDITIONER_H
#define PARCSR_ILUAMG_PRECONDITIONER_H

class ParCSR_ILUAMG_Preconditioner : public ParCSR_Preconditioner
{
private:
	int num_levels;
	ParCSR_Matrix* A;
	ParCSR_Matrix* P;
	ParCSR_Matrix* R;
	class Smoother;
	Smoother* S;

	Par_Vector* r;
	Par_Vector* e;
	Par_Vector* g;

	void Strength(const ParCSR_Matrix& A, bool* LocLoc_s, bool* LocExt_s) const;
	void Coarsening(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map, int AGGRESSIVE = 0) const;
	void Interpolation(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map, ParCSR_Matrix& P, int AGGRESSIVE = 0) const;
	void RAP(const ParCSR_Matrix& R, const ParCSR_Matrix& A, const ParCSR_Matrix& P, ParCSR_Matrix& C) const;
	void Sparsification(const ParCSR_Matrix& A, ParCSR_Matrix& B) const;
	void Smooth(int level, const Par_Vector& b, const Par_Vector& x) const;
	void Restriction(int level, const Par_Vector& b, const Par_Vector& x, const Par_Vector& g) const;
	void Prolongation(int level, const Par_Vector& e, const Par_Vector& x) const;
	void V_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;
	void W_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;
	void F_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;

public:
	int    MaxLevels;
	int    CoarseSize;
	double StrengthThreshold;
	int    AggNumLevels;
	int    CoarsenType;
	int    InterpMinElements;
	int    InterpMaxElements;
	double TruncationFactor;
	double SparsificationThreshold;
	int    CycleType;
	int    PreSweeps;
	int    PostSweeps;
	int    CoarseSweeps;
	int    PrintStats;

	ParCSR_ILUAMG_Preconditioner();
	ParCSR_ILUAMG_Preconditioner(
	int    max_levels,
	int    coarse_size,
	double strength_threshold,
	int    agg_num_levels,
	int    coarsen_type,
	int    interp_min_elements,
	int    interp_max_elements,	
	double truncation_factor,
	double sparsification_threshold,
	int    cycle_type,
	int    pre_sweeps,
	int    post_sweeps,
	int    coarse_sweeps,
	int    print_stats);
	~ParCSR_ILUAMG_Preconditioner();

	void Free();
	void Setup(const ParCSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef BSR_MATRIX_H
#define BSR_MATRIX_H

struct BSR_Matrix : public Operator
{
	bool ref;
	int size[2];
	int block_size;

	int* rowptr;
	int* colind;
	double* values;

	BSR_Matrix();
	BSR_Matrix(int n, int m, int bsize, int* rowptr, int* colind, double* values, int ref);
	BSR_Matrix(const BSR_Matrix& A);
	BSR_Matrix(BSR_Matrix&& A);
	~BSR_Matrix();
	BSR_Matrix& operator=(const BSR_Matrix& A);
	BSR_Matrix& operator=(BSR_Matrix&& A);

	void Free();
	void Reference(const BSR_Matrix& A);
	int InSize() const;
	int OutSize() const;
	void Apply(const Vector& x, const Vector& y) const;
};

#endif

#ifndef PARBSR_MATRIX_H
#define PARBSR_MATRIX_H

struct ParBSR_Matrix : public Par_Operator
{
	BSR_Matrix LocLoc;
	BSR_Matrix LocExt;

	int num_neighbors;
	int* neighbor_ranks;

	int* Ext_ptr;
	int* Ext_map;

	int* send_ptr;
	int* send_list;
	double* send_buffer;
	
	Vector x_Ext;

	ParBSR_Matrix();
	ParBSR_Matrix(MPI_Comm comm, int LocLoc_rows, int LocLoc_cols, int LocExt_cols, int block_size, 
	int* LocLoc_rowptr, int* LocLoc_colind, double* LocLoc_values, int LocLoc_ref, 
	int* LocExt_rowptr, int* LocExt_colind, double* LocExt_values, int LocExt_ref,
	int num_neighbors, int* neighbor_ranks, int* Ext_ptr, int* Ext_map);
	ParBSR_Matrix(const ParBSR_Matrix& A);
	ParBSR_Matrix(ParBSR_Matrix&& A);
	~ParBSR_Matrix();
	ParBSR_Matrix& operator=(const ParBSR_Matrix& A);
	ParBSR_Matrix& operator=(ParBSR_Matrix&& A);
	
	void Free();
	void Reference(const ParBSR_Matrix& A);
	void ExchangeHalo(const Vector& x_Loc, const Vector& x_Ext) const;
	void SetupHalo();
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& x, const Par_Vector& y) const;
};

#endif

#ifndef PARBSR_MATVEC_H
#define PARBSR_MATVEC_H

void ParBSR_MatVec(double alpha, const ParBSR_Matrix& A, const Par_Vector& x, double beta, const Par_Vector& y, const Par_Vector& w);
void ParBSR_MatVec(const ParBSR_Matrix& A, const Par_Vector& x, const Par_Vector& y);

#endif

#ifndef PARBSR_GENERATOR_H
#define PARBSR_GENERATOR_H

class ParBSR_Generator
{
public:
	MPI_Comm comm;

	ParBSR_Generator();
	ParBSR_Generator(MPI_Comm comm);
	virtual ~ParBSR_Generator();
	virtual void operator()(ParBSR_Matrix& A) const = 0;
};

#endif

#ifndef PARBSR_LAPLACIAN_GENERATOR_H
#define PARBSR_LAPLACIAN_GENERATOR_H

class ParBSR_Laplacian_Generator : public ParBSR_Generator
{
public:
	int nx;
	int ny;
	int nz;

	int bx;
	int by;
	int bz;

	int Px;
	int Py;
	int Pz;

	double cx;
	double cy;
	double cz;

	ParBSR_Laplacian_Generator();
	ParBSR_Laplacian_Generator(MPI_Comm comm);
	ParBSR_Laplacian_Generator(MPI_Comm comm, int nx, int ny, int nz, int bx, int by, int bz, int Px, int Py, int Pz, double cx, double cy, double cz);
	void operator()(ParBSR_Matrix& A) const;
};

#endif

#ifndef PARBSR_LAPLACIAN27PT_Generator_H
#define PARBSR_LAPLACIAN27PT_Generator_H

class ParBSR_Laplacian27pt_Generator : public ParBSR_Generator
{
public:
	int nx;
	int ny;
	int nz;

	int bx;
	int by;
	int bz;

	int Px;
	int Py;
	int Pz;

	ParBSR_Laplacian27pt_Generator();
	ParBSR_Laplacian27pt_Generator(MPI_Comm comm);
	ParBSR_Laplacian27pt_Generator(MPI_Comm comm, int nx, int ny, int nz, int bx, int by, int bz, int Px, int Py, int Pz);
	void operator()(ParBSR_Matrix& A) const;
};

#endif

#ifndef PARBSR_IJ_GENERATOR_H
#define PARBSR_IJ_GENERATOR_H

class ParBSR_IJ_Generator : public ParBSR_Generator
{
public:
	int Loc_rows;
	int block_size;
	const long* col_partitioning;
	const int* LocGlo_rowptr;
	const long* LocGlo_colind;
	const double* LocGlo_values;
	int pattern_symmetry;

	ParBSR_IJ_Generator();
	ParBSR_IJ_Generator(MPI_Comm comm);
	ParBSR_IJ_Generator(MPI_Comm comm, int Loc_rows, int block_size, const long* col_partitioning, const int* LocGlo_rowptr, const long* LocGlo_colind, const double* LocGlo_values, int pattern_symmetry);
	void operator()(ParBSR_Matrix& A) const;
};

#endif

#ifndef PARBSR_PRECONDITIONER_H
#define PARBSR_PRECONDITIONER_H

class ParBSR_Preconditioner : public Par_Operator
{
public:
	virtual ~ParBSR_Preconditioner();
	virtual void Setup(const ParBSR_Matrix& A, int REUSE = 0) = 0;
};

#endif

#ifndef PARBSR_JACOBI_PRECONDITIONER_H
#define PARBSR_JACOBI_PRECONDITIONER_H

class ParBSR_Jacobi_Preconditioner : public ParBSR_Preconditioner
{
private:
	int size;
	int block_size;
	double* D_LU;

public:
	ParBSR_Jacobi_Preconditioner();
	~ParBSR_Jacobi_Preconditioner();
	void Free();
	void Setup(const ParBSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARBSR_SOR_PRECONDITIONER_H
#define PARBSR_SOR_PRECONDITIONER_H

class ParBSR_SOR_Preconditioner : public ParBSR_Preconditioner
{
private:
	double* D_Loc;
	double* D_LU_Loc;
	BSR_Matrix L_LocLoc;
	BSR_Matrix L_LocExt;
	BSR_Matrix U_LocLoc;
	BSR_Matrix U_LocExt;
	Vector x_Ext;

	bool* interior;

	int label;
	int num_neighbors;
	int* neighbor_labels;
	int* neighbor_ranks;

	int* recv_ptr;
	int* send_ptr;
	int* send_list;
	double* send_buffer;

public:
	int Type;
	double Factor;

	ParBSR_SOR_Preconditioner();
	ParBSR_SOR_Preconditioner(int type, double factor);
	~ParBSR_SOR_Preconditioner();

	void Free();
	void Setup(const ParBSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARBSR_ILU_PRECONDITIONER_H
#define PARBSR_ILU_PRECONDITIONER_H

class ParBSR_ILU_Preconditioner : public ParBSR_Preconditioner
{
private:
	double* D_LU_Loc;
	BSR_Matrix L_LocLoc;
	BSR_Matrix L_LocExt;
	BSR_Matrix U_LocLoc;
	BSR_Matrix U_LocExt;
	Vector x_Ext;

	bool* interior;

	int label;
	int num_neighbors;
	int* neighbor_labels;
	int* neighbor_ranks;

	int* recv_ptr;
	int* send_ptr;
	int* send_list;
	double* send_buffer;

public:
	ParBSR_ILU_Preconditioner();
	~ParBSR_ILU_Preconditioner();
	
	void Free();
	void Setup(const ParBSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARBSR_BLOCKJACOBISOR_PRECONDITIONER_H
#define PARBSR_BLOCKJACOBISOR_PRECONDITIONER_H

class ParBSR_BlockJacobiSOR_Preconditioner : public ParBSR_Preconditioner
{
private:
	int num_threads;
	double* D;
	double* D_LU;
	BSR_Matrix L;
	BSR_Matrix U;

public:
	int Type;
	double Factor;

	ParBSR_BlockJacobiSOR_Preconditioner();
	ParBSR_BlockJacobiSOR_Preconditioner(int type, double factor);
	~ParBSR_BlockJacobiSOR_Preconditioner();

	void Free();
	void Setup(const ParBSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARBSR_BLOCKJACOBIILU_PRECONDITIONER_H
#define PARBSR_BLOCKJACOBIILU_PRECONDITIONER_H

class ParBSR_BlockJacobiILU_Preconditioner : public ParBSR_Preconditioner
{
private:
	int num_threads;
	double* D_LU;
	BSR_Matrix L;
	BSR_Matrix U;

public:
	ParBSR_BlockJacobiILU_Preconditioner();
	~ParBSR_BlockJacobiILU_Preconditioner();
	void Free();
	void Setup(const ParBSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARBSR_AMG_PRECONDITIONER_H
#define PARBSR_AMG_PRECONDITIONER_H

class ParBSR_AMG_Preconditioner : public ParBSR_Preconditioner
{
private:
	int num_levels;
	ParBSR_Matrix* A;
	ParCSR_Matrix* P;
	ParCSR_Matrix* R;
	class Smoother;
	Smoother* S;

	Par_Vector* r;
	Par_Vector* e;
	Par_Vector* g;

	void Strength(const ParCSR_Matrix& A, bool* LocLoc_s, bool* LocExt_s) const;
	void Coarsening(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map, int AGGRESSIVE = 0) const;
	void Interpolation(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map, ParCSR_Matrix& P, int AGGRESSIVE = 0) const;
	void RAP(const ParCSR_Matrix& R, const ParBSR_Matrix& A, const ParCSR_Matrix& P, ParBSR_Matrix& C) const;
	void Sparsification(const ParBSR_Matrix& A, ParBSR_Matrix& B) const;
	void Smooth(int level, const Par_Vector& b, const Par_Vector& x) const;
	void Restriction(int level, const Par_Vector& b, const Par_Vector& x, const Par_Vector& g) const;
	void Prolongation(int level, const Par_Vector& e, const Par_Vector& x) const;
	void V_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;
	void W_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;
	void F_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;

public:
	int    MaxLevels;
	int    CoarseSize;
	double StrengthThreshold;
	int    AggNumLevels;
	int    CoarsenType;
	int    InterpMinElements;
	int    InterpMaxElements;
	double TruncationFactor;
	double SparsificationThreshold;
	int    CycleType;
	int    PreSweeps;
	int    PostSweeps;
	int    CoarseSweeps;
	int    RelaxationType;
	double RelaxationFactor;
	int    PrintStats;

	ParBSR_AMG_Preconditioner();
	ParBSR_AMG_Preconditioner(
	int    max_levels,
	int    coarse_size,
	double strength_threshold,
	int    agg_num_levels,
	int    coarsen_type,
	int    interp_min_elements,
	int    interp_max_elements,	
	double truncation_factor,
	double sparsification_threshold,
	int    cycle_type,
	int    pre_sweeps,
	int    post_sweeps,
	int    coarse_sweeps,
	int    relaxation_type,
	double relaxation_factor,
	int    print_stats);
	~ParBSR_AMG_Preconditioner();
	
	void Free();
	void Setup(const ParBSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef PARBSR_ILUAMG_PRECONDITIONER_H
#define PARBSR_ILUAMG_PRECONDITIONER_H

class ParBSR_ILUAMG_Preconditioner : public ParBSR_Preconditioner
{
private:
	int num_levels;
	ParBSR_Matrix* A;
	ParCSR_Matrix* P;
	ParCSR_Matrix* R;
	class Smoother;
	Smoother* S;

	Par_Vector* r;
	Par_Vector* e;
	Par_Vector* g;

	void Strength(const ParCSR_Matrix& A, bool* LocLoc_s, bool* LocExt_s) const;
	void Coarsening(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map, int AGGRESSIVE = 0) const;
	void Interpolation(const ParCSR_Matrix& A, const bool* LocLoc_s, const bool* LocExt_s, int* CF_map, ParCSR_Matrix& P, int AGGRESSIVE = 0) const;
	void RAP(const ParCSR_Matrix& R, const ParBSR_Matrix& A, const ParCSR_Matrix& P, ParBSR_Matrix& C) const;
	void Sparsification(const ParBSR_Matrix& A, ParBSR_Matrix& B) const;
	void Smooth(int level, const Par_Vector& b, const Par_Vector& x) const;
	void Restriction(int level, const Par_Vector& b, const Par_Vector& x, const Par_Vector& g) const;
	void Prolongation(int level, const Par_Vector& e, const Par_Vector& x) const;
	void V_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;
	void W_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;
	void F_Cycle(int level, const Par_Vector& b, const Par_Vector& x) const;

public:
	int    MaxLevels;
	int    CoarseSize;
	double StrengthThreshold;
	int    AggNumLevels;
	int    CoarsenType;
	int    InterpMinElements;
	int    InterpMaxElements;
	double TruncationFactor;
	double SparsificationThreshold;
	int    CycleType;
	int    PreSweeps;
	int    PostSweeps;
	int    CoarseSweeps;
	int    PrintStats;

	ParBSR_ILUAMG_Preconditioner();
	ParBSR_ILUAMG_Preconditioner(
	int    max_levels,
	int    coarse_size,
	double strength_threshold,
	int    agg_num_levels,
	int    coarsen_type,
	int    interp_min_elements,
	int    interp_max_elements,	
	double truncation_factor,
	double sparsification_threshold,
	int    cycle_type,
	int    pre_sweeps,
	int    post_sweeps,
	int    coarse_sweeps,
	int    print_stats);
	~ParBSR_ILUAMG_Preconditioner();
	
	void Free();
	void Setup(const ParBSR_Matrix& A, int REUSE = 0);
	int InSize() const;
	int OutSize() const;
	void Apply(const Par_Vector& b, const Par_Vector& x) const;
};

#endif

#ifndef YHAMG_H
#define YHAMG_H

typedef Vector                               Vec;
typedef MultiVector                          MultiVec;
typedef Par_Vector                           ParVec;
typedef Par_MultiVector                      ParMultiVec;

typedef Operator                             OpBase;
typedef Par_Operator                         ParOpBase;
typedef Par_Zero_Operator                    ParOpZero;
typedef Par_Identity_Operator                ParOpIdent;
typedef Par_Compound_Operator                ParOpComp;
typedef Par_Chebyshev_Polynomial             ParOpChebyshev;

typedef Par_Solver                           ParSolverBase;
typedef Par_CG_Solver                        ParSolverCG;
typedef Par_GMRES_Solver                     ParSolverGMRES;
typedef Par_BiCGStab_Solver                  ParSolverBCGS;
typedef Par_PipeCG_Solver                    ParSolverPCG;
typedef Par_PipeGMRES_Solver                 ParSolverPGMRES;
typedef Par_PipeBiCGStab_Solver              ParSolverPBCGS;

typedef Par_EigenSolver                      ParEigenSolverBase;
typedef Par_LOBPCG_EigenSolver               ParEigenSolverLOBPCG;

typedef COO_Matrix                           COO;
typedef CSR_Matrix                           CSR;
typedef ParCSR_Matrix                        ParCSR;

typedef ParCSR_Preconditioner                ParCSRPrecondBase;
typedef ParCSR_Jacobi_Preconditioner         ParCSRPrecondJacobi;
typedef ParCSR_SOR_Preconditioner            ParCSRPrecondSOR;
typedef ParCSR_ILU_Preconditioner            ParCSRPrecondILU;
typedef ParCSR_BlockJacobiSOR_Preconditioner ParCSRPrecondBJSOR;
typedef ParCSR_BlockJacobiILU_Preconditioner ParCSRPrecondBJILU;
typedef ParCSR_AMG_Preconditioner            ParCSRPrecondAMG;
typedef ParCSR_ILUAMG_Preconditioner         ParCSRPrecondILUAMG;

typedef ParCSR_Generator                     ParCSRGeneratorBase;
typedef ParCSR_IJ_Generator                  ParCSRGeneratorIJ;
typedef ParCSR_Laplacian_Generator           ParCSRGeneratorLaplace;
typedef ParCSR_Laplacian27pt_Generator       ParCSRGenerator27pt;

typedef BSR_Matrix                           BSR;
typedef ParBSR_Matrix                        ParBSR;

typedef ParBSR_Preconditioner                ParBSRPrecondBase;
typedef ParBSR_Jacobi_Preconditioner         ParBSRPrecondJacobi;
typedef ParBSR_SOR_Preconditioner            ParBSRPrecondSOR;
typedef ParBSR_ILU_Preconditioner            ParBSRPrecondILU;
typedef ParBSR_BlockJacobiSOR_Preconditioner ParBSRPrecondBJSOR;
typedef ParBSR_BlockJacobiILU_Preconditioner ParBSRPrecondBJILU;
typedef ParBSR_AMG_Preconditioner            ParBSRPrecondAMG;
typedef ParBSR_ILUAMG_Preconditioner         ParBSRPrecondILUAMG;

typedef ParBSR_Generator                     ParBSRGeneratorBase;
typedef ParBSR_IJ_Generator                  ParBSRGeneratorIJ;
typedef ParBSR_Laplacian_Generator           ParBSRGeneratorLaplace;
typedef ParBSR_Laplacian27pt_Generator       ParBSRGenerator27pt;

#endif

}